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NRQCD results for Upsilon and Charmonium using both dynamical and quenched configurations are presented. 
We investigate dependence on the light dynamical quark mass. Preliminary dynamical (n/ ~ 2) Charmonium 
data are combined with quenched results to extract the strong coupling constant for the physical number 

of light dynamical quarks, Uf — 3. Good agreement is found with calculations based on the Upsilon system. We 
show that a discrepancy in found between the Upsilon and Charmonium systems in the quenched theory, 

disappears upon extrapolating to the physical number of flavors. Results for the strong coupling constant a-jgj 
are presented and sources of systematic error investigated. 



1. Introduction 

The Nonrelativistic QCD (NRQCD) ap- 
proach to heavy quarks on the lattice has been 
applied very successfully to the T and Charmo- 
nium systems in recent years Good agree- 
ment is found between lattice simulations and ob- 
served quarkonium levels. These investigations 
have also enabled us to determine the b-quark 
pole- and MS'-masses Q and the strength of the 
strong coupling constant, 0. 

During the past year we have continued studies 
of quarkonium systems with particular emphasis 
on refining our as determination. We have : 

1. a new high statistics T study on 16'^ x 32 
Kogut et al configurations with = 
at /3 = 6.0. 

2. increased statistics for T studies on 16^ x 32 
HEMCGC g configurations with n/ = 2, 
arriq = 0.01 staggered dynamical quarks at 
P = 5.6. 

3. new T results using HEMCGC dynamical 
configurations with aruq = 0.025. 



4. new ( preliminary) Charmonium data using 
16^ X 24 dynamical MILC [|j configurations 
(n/ = 2, arriq = 0.0125) at /3 = 5.145. 

This has led to the following main results , 

• an T spectroscopy update 

• a look at the mg-dependence of a^^"^'' 

• extrapolation of ap\ ap' a'p^~'^^ from 
both T and Charmonium. 

a 5 ~ 6 sigma discrepancy at = is 
removed through unquenching and extrap- 
olation to the physical = 3. 

• establishment of a reliable calculation for 
^M5 mass. 

2. T in the Quenched Approximation 

The NRQCD action, details of how heavy 
quark propagators are evaluated, the smearing 
functions and the fitting procedures used are de- 
scribed in detail in several publications 
We will review them briefly. 
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With ~ 0.1, the b quarks in the T sys- 
tem are quite nonrelativistic. The spUttings be- 
tween spin-averaged levels are around 500GeV 
[©(Mfet)^)], which is much smaller than the T 
mass [0(2M{,)], indicating that a systemic expan- 
sion of the QCD Hamiltonian in powers of is 
appropriate. The continuum action density, cor- 
rect through ©(Mfc^;'*), is broken down according 
to 



He 



fl) 



ifg°"* and SH'^""'* are given explicitly in ref. , we 
give here their lattice counterparts Hq and SH: 



with 

Gi(y)^^Gi(y-x)r(-)(x)e'^^-'' 

X 

where the trace is over spin and color. r'^*'^\ the 
smearing at the source is distinguished from that 
at the sink, r('*'=) . 

Two main fitting procedures were used to ex- 
tract energies and amplitudes. For the first pro- 
cedure, a matrix of correlations with riscj^sfc = 
1,2,3 for S states and nsc,nsk = 1, 2 for spin sin- 
glet P states was fit simultaneously, each correla- 
tion being fit to an ansatz of the form 
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and 
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The last two terms in SH come from finite lattice 
spacing corrections to the lattice Laplacian and 
lattice time derivative, respectively. 

The quark propagators are determined from 
evolution equations that specify the propagator 
value, for t > 0, in terms of the value on the pre- 
vious timeslice; 



G 



t+i = 



1 



2n 



ul\ 1 



2n 



(1 - aSH)Gt 



The quark propagators are combined with the 
smearing operator F to produce a meson propa- 
gator 



Gqq{nsc,nsk;t) = ^ a{nsc,k)a*{nsk,k)c 



-E,.,t 



fe=i 



The second procedure involved fitting a row of 
smeared-local correlations, Hsc — 1,2,3 or rise = 
1,2 with Usk — loc, simultaneously, each correla- 
tion being fit to an ansatz of the form 



fe=l 



A comparison of the dimensionless IS-l-^P and 
1S-2S T splittings, once C(a^) gluonic corrections 
are made (using perturbation theory aSEg ~ 
0.0036 for the IS and ~ 0.0023 for the 2S 
level), with experimental values for these split- 
tings yields the following values for a^^ 



IS - IP: 
IS - 2S: 



2.59(6) GeV 
2.34(5) GeV 



These results represent an increased statistical ac- 
curacy of a factor of 2 over our previous results 
Taking an average of a^^ = 2.4 GeV, the 
updated NRQCD T spectrum is shown in figures 
1 & 2 (see section 3). The agreement between 
figures 1 & 2 and their previous counterparts Q 
is good - a significant shift {2a) is only observed 
for the 3S state. This may be an indication that 
previous statistical errors, which included an ele- 
ment of fitting uncertainty, were overestimates. 
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3. T with Dynamical Quarks 



We now have twice the statistics on nt = 2 
configurations as compared to a year ago and 
have accumulated 6,400 meson propagators per 
channel and per choice of smearing at source and 
sink. As in the past we find it crucial to use 
multi-exponential fits to several correlations si- 
multaneously, in order to extract ground state 
and one or two excited state energies for the S- 
states and the spin averaged P-states. We re- 
fer to |^,|^ for details and mention here just that 
with higher statistics we are now much more sen- 
sitive to systematics. To give one example, the 
statistical errors of the dimensionless IS-IP split- 
ting are of order ~ 0.005 which is comparable to 
shifts expected from O(a^) corrections to the glu- 
onic action. The latter have been calculated per- 
turbatively to be aSEg ^ 0.0057 for the IS and 
~ 0.0034 for the 2S level. With several groups 
now working with improved gluonic actions 
one would eventually like to check aSEg nonper- 
turbatively. These corrections are crucial when 
we check for scaling by comparing T simulation 
results at different P values. 

By comparing the dimensionless IS-IP and IS- 
2S T splittings with experiment, one determines 



GeV 
10.5- 



We find 

IS - IP: 
IS - 2S: 



= 2.44(7) GeV 
= 2.37(10) GeV 



Using an average a~^=2.4GeV we show our up- 
dated results for the NRQCD T spectrum in fig- 
ures 1 & 2. Upon focusing on the IP and 2S 
levels in figure 1, one notices better agreement 
with experiment in the dynamical theory than in 
the quenched theory. This reflects the fact that 
a~^'s from the IS-IP and 1S-2S splittings dis- 
agree with each other (at the ^ 4-sigma level) in 
the quenched theory. As discussed in reference 
[t] the effect of this discrepancy in disappears 
when one extrapolates to the ny = 3 theory. This 
remains true also of our new nj = and Uf ~ 2 
data. A similar phenomenon of different quan- 
tities giving different a^^'s in the quenched the- 
ory, whose effects then disappear for the physi- 
cal number of dynamical flavors, will be discussed 
and explained in section 5 where we compare T 
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Figure 1. T Spectrum for nj — (open circles) 
and Uf = 2 (full circles). In both cases, a^^ = 
2.4GeV was used. The dashed horizontal lines 
denote experimental values. 



and Charnionium results. 

Figure 2 shows spin splittings. One sees that 
there are some differences between experiment 
and simulations in the Xb splittings. This is par- 
ticularly noticeable in the quenched case and one 
is tempted to try to extrapolate in n/. Before 
carrying out that exercise, we need to reexamine 
our fine structure fits. The fitting procedure used 
for the P-state fine structure differs from what 
was done for the S-states and the ^Pi levels to 
get the IS-IP and 1S-2S spfittings. For the latter 
levels we used multi-exponential multi-correlation 
fits as mentioned above. In order to extract the 
small fine-structure splittings (these are of order 
10 MeV = size of the statistical errors in the ^Pi 
level) we looked at ratios between "^Pj and ^Pi 
correlations. These were fit to a single exponen- 
tial. In past work it was not possible to find a 
signal for excited state contributions to these ra- 
tios. With our current higher statistics data it 
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Figure 2. T SPIN SPLITTINGS : Symbols 
have the same meanings as in Figure 1. 



may be possible and necessary to take such con- 
tributions into account. Further work is required 
before one can sort out quenching versus relativis- 
tic corrections to the fine structure splittings. 

The last topic to discuss in this section on 
dynamical T data, concerns dependence on the 
mass of the dynamical quarks. The data dis- 
cussed above used HEMCGC dynamical config- 
urations with Uf — 2 staggered fermions of di- 
mensionless mass amq = 0.01. With inverse lat- 
tice spacing ~ 2.4GeV, this corresponds to 
light quarks considerably heavier than the up- or 
down- quarks. This should not pose a problem as 
long as niq << typical momenta in heavy-heavy 
systems of about 0.5 ^ IGeV (with the higher 
number applicable for T). Hence heavy quark 
physics with dynamical light quarks should be 
much less sensitive to extrapolations in the light 
quark mass than light quark physics. We have 
made the first attempt to test this hypothesis by 
repeating our T runs with aruq — 0.025 dynami- 
cal HEMCGC configurations. The inverse lattice 
spacings are found to be 



anig = 0.025 



IS 
IS 



IP: 
2S: 



fl-i = 2.28(15) GeV 
= 2.17(12) GeV 



configurations, one finds a 1-sigma difference for 
the S-P and a 2-sigma change for the 1S-2S a~^. 
It is not clear with present statistics, whether we 
are observing a true systematic effect, although 
the shift downwards in the central values would 
be in the right direction for an aniq that has be- 
come too large. At given /3 increasing the dynam- 



ical 



will reduce the difference in /3 between 



that and a quenched simulation of the same a~^. 

An accurate estimate of is the crucial input 
for lattice determinations of the strong coupling 
constant ||^,^ . So what one is really interested in 
here is the dependence of Us extracted from 
lattice quarkonium studies. We use an ap defined 
through the plaquette value, 



47r 3 41 

-InWi.i = yap(— — )[l-(1.185 + 0.070n/)ap] 

with a^^ from quarkonium splittings setting the 
scale. Using the above a~^'s and evolving the 
couplings perturbatively to a common reference 
scale of 8.2GeV, one finds 



S-P 



„(»/=2)r 



a 



("/=2)r 



1S-2S 



a. 



a 



(n/=2)r 



5.2GeV] = 0.1793(16) am, = 0.01 
5.2GeV] = 0.1760(35) am, = 0.025 



5.2GeV] = 0.1777(23) am, = 0.01 
5.2GeV] = 0.1735(28) am, = 0.025 



One could try extrapolating to m, — > 0. Pertur- 
bation theory tells us to extrapolate quadratically 
in TO,. One then finds results very close to the 
arUq = 0.01 values (0.1799(20) and 0.1785(28) re- 
spectively from S-P and 1S-2S) More studies of 
the m,-dependence in dynamical simulations of 
quarkonium systems are called for, with higher 
statistics and several to, values. At the moment 
our calculations indicate that the effects of ex- 
trapolating from am, = 0.01 down to realistic 
light quark masses are very small and less than 
our statistical errors. Hence, we will continue to 
use the and ap values from the am, = 0.01 
simulations for which we have the best statistics. 



Comparing with the a 



s from the am, = 0.01 
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4. Us Determination and Investigation of 
Systematic Errors 

One of the simplest quantities to calculate in 
lattice QCD, and to yield a value for the strong 
coupling constant, is the expectation value of the 
1x1 Wilson loop operator. The coupling con- 
stant, ap, is defined through the plaquette value 
by the perturbative relation in section 3. With 
this definition, ap coincides through 0{{ap)^) 
with the coupling ay, defined in refs. pO| , pl| . 

Using the values calculated in sections 2 & 
3 [aniq = 0.01 results from section 3) to set the 
scale and evolving the couplings perturbatively 
using the 2-loop formula to the common reference 
scale of 8.2 GeV, one obtains 



IS - IP 

„,("/=0)r 



a 



(n/=2) 



5.2GeV] = 0.1551(11) 
[8.2GeV] = 0.1793(16) 



IS - 2S 

X. p 



'[8.2GeV] = 0.1505(9) 



a^"^-^^[8.2GeV] = 0.1777(23) 

To obtain physical results, one must extrapolate 
the couplings to ti/ = 3. The inverse of the cou- 
pling is known to be almost linear for 
small changes in n/, hence extrapolating the in- 
verse couplings one obtains 



a 



5.2GeV] = 0.1945(30) IS-IP 
5.2GeV] = 0.1953(43) 1S-2S 



The concordance of these values is more readily 
appreciated from figure 5. 

To enable comparison with other determina- 
tions we convert to the MS scheme through the 
relation 



MS 



') ) 



(2) 



The factor of exp(5/6) absorbs the n/ depen- 
dence in (at least) the term. Substituting for 
the a'p'~^^ values shown above, evolving pertur- 



batively to a scale of 1.3GeV, matching from 3 
to 4 flavors of dynamical quarks, evolving once 
more to 4.3GeV, matching from 4 to 5 flavors 
of dynamical quarks and finally evolving to the 
Mzo = 91.2GeV scale, we obtain 



■msWzA = 0.1152(24) IS-IP 
-MsWzA = 0.1154(26) 1S-2S 



The error in these numbers is dominated by the 
terms in equation 2. The central value of 
course has zero for this coefficient; one standard 
deviation allows this coefficient to have magni- 
tude 1. Ref. ||l^ has shown that for an = 

theory, the coefficient of the (a^^') term in 
equation ^ is 0.96, giving a central value for 
a-^[Mzo] of 0.117. The remaining uncertainty 
would then be the n/ dependence of the co- 
efficient, and that may be considerably less than 
the previous uncertainty. 

Although this coefficient is the main source of 
error, we have also tested the robustness of our 
result to other factors. In one such test, the 0{a?') 
gluonic corrections to the action were omitted, 
yielding the results 

ajjg[Mz«] = 0.1145(24) IS-IP 
a-^[Mz«] = 0.1154(26) 1S-2S 

The 1S-2S value is identical to that obtained pre- 
viously with the corrections in place since the cor- 
rections to the IS and 2S levels partially cancel 
each other. The IS- IP value is slightly lower be- 
cause there are no P-state 0{a?) gluonic correc- 
tions. So, as expected, although the spectrum re- 
sults are accurate enough to distinguish the pres- 
ence of this correction our final value for a is not. 

The effect of tadpole-improvement was tested 
by using quenched results from with ua — \. 
We have no dynamical results without tadpole- 
improvement so we modified our existing dynami- 
cal results to simulate this situation. We adjusted 
them by the difference of the tadpole-improved 
and tadpole-unimproved quenched energies, since 
the effect of unquenching on this difference is a 
second order effect. The subsequent values 
from the 1S-2S splitting gave rise to the result 

ajjg[Mz«] = 0.1149(27) 1S-2S 
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Figure 3. Variation of cejfg with threshold at 
which matching from 3 to 4 dynamical quark fla- 
vors takes place. 



Figure 4. Variation of cejfg with threshold at 
which matching from 4 to 5 dynamical quark fla- 
vors takes place. 



This result is very close to that obtained with 
tadpole-improvement since this has little effect on 
spin-independent splittings (it is crucial for spin 
splittings). 

The effect of the correction terms, SH, in the 
action of equation^ was tested using quenched S- 
state energies from with 6H — and the dy- 
namical S-state energies used previously but with 
a modification equivalent to that described above 
for the omission of tadpole-improvement. This 
time we obtained the result 

ajjs[Mzo] = 0.1165(27) 1S-2S 

Although the value for cxjjg is larger than that 
for non-zero SH, the increase is masked by the 
uncertainty incurred when we converted from ap 
to ajj-g. It is clear that higher order relativis- 
tic corrections that we do not include would be 
completely invisible. 

The final factors we investigated were the 
thresholds at which we matched from 3 to 4 and 
from 4 to 5 flavors of dynamical quarks . Fig- 
ure 3 shows values of ctjjg obtained from the IS- 



IP splitting where the lower matching threshold 
(3 to 4 flavors) was varied around the l.SGeV 
value while the upper matching threshold (4 to 
5 flavors) was fixed at 4.3GeV. Figure 4 is simi- 
lar; this time the upper matching threshold was 
varied around 4.3GeV with the lower fixed at 
l.SGeV. For simplicity we used the values from 
ref. |l| of 1.3 GeV and 4.3 GeV for the quark 
masses, mq{niq). As can be seen from the figures 
there are comfortable ranges of at least 1.5GeV 
where there is essentially no effect on our final 
result, in accord with the findings of ref. [ p2| . 

5. Charmonium with Dynamical Quarks 
and a New Estimate of as 

The that feeds into lattice determinations 
of as can, in principle, be obtained from a vari- 
ety of physical quantities. The reason level split- 
tings in the T system are optimal, is that sta- 
tistical errors can be reduced rather efficiently 
in these systems, and because systematic errors 
are also under good control. Another similar sys- 
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tern that can be used for as determinations is the 
Charmonium system Our collaboration has 
studied Charmonium in the quenched approxima- 
tion using NRQCD heavy fermions on /3 = 5.7 
UKQCD configurations. These results are cov- 
ered by Christine Davies at this conference 
We now also have preliminary dynamical Char- 
monium data on 71/ = 2 MILC ^ configurations 
at /3 = 5.145. Using the Charmonium spin aver- 
aged S-P splitting to determine , one ends up 



with a value for a 



(nf=2) 



of, 



a 



(»/=2) 



[8.2GeV] = 0.1758(36) Charm 



which should be compared with the T value of 
0.1793(16) from the previous section. 

ttp need not agree between Charmonium and T 
for 71/ = and Uf — 2. However there should be 
agreement for the physical number of dynamical 
fiavors, 71/ = 3. In fact based on the quenched re- 
sults in section 4 and those discussed by C. Davies 
jll one finds for ri/ = 



a^"^-°^[8.2GeV] = 0.1480(13) 
a^"^=°^[8.2GeV] =0.1551(11) 



Charm 
T 



We believe the 5 ~ 6 sigma discrepancy between 
Charmonium and T is largely a quenching effect 
combined with the fact that the two quarkonium 
systems have different characteristic energy scales 
and (g^ > Qc)- ^ quenched theory a{q) 
will run incorrectly between the two g*'s. As a 
consequence the T S-P splitting will be under- 
estimated relative to the same splitting in Char- 
monium. Upon comparing with real experimen- 
tal data the T simulation results will lead to a 
larger than in Charmonium and hence also 
to a larger ap[8.2GeU]. 

The 71/ = 2 a^^~^'''s agree better between 
Charmonium and T, although one does not ex- 
pect complete agreement until one has extrapo- 
lated to the physical number of flavors ri/ = 3. 
We extrapolate in 1/ap and obtain, 



a^"^-^^[8.2GeV] = 0.1940(67) 
a^^"^'[8.2GeV] = 0.1945(30) 



Charm 
T 



So, extrapolation to the correct number of flavors 
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Figure 5. ap^ [8.2GeV] versus Nf with the 
scale set by the T IS-IP (full chcles), the T 1S-2S 
(open circles) or the Charmonium IS-IP (boxes) 
splitting. 



has removed the discrepancy altogether. The ii/ 
dependence is summarized in figure 5. The errors 
on the Charmonium results are still large, but 
we are working towards improving them. The 
^(nj,_3) Charmonium and the updated value 
from T agree well with our previous published 
value 0. 

In figure 5 we also show ap^'' from T 1S-2S 
splittings (the open circles) . Again a discrepancy 
at 71/ = disappears upon going to the ii/ = 3 
theory. Arguing why, in the same quarkonium 



system, one expects a larger a^^""^ from a IS- 
IP determination than from the 1S-2S splitting, 
is somewhat subtle. It no longer suffices to talk 
about just one typical characteristic as we did 
when comparing T with Charmonium. One must 
look more carefully into the details of individual 
binding energies for the S- and P-states and the 
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a{q) values associated with them. As discussed 
in reference 0, in a theory with incorrect run- 
ning of a (a running that is too rapid) the IS-IP 
spUtting will again be underestimated relative to 
the 1S-2S splitting, leading to the phenomenon 
observed in figure 5. The important point here 
is that we have now observed subtle quenching 
effects in several ways and that in each case such 
effects are removed by including the appropriate 
number of dynamical light quarks. 

6. Summary 

We have carried out several tests/updates of 
lattice determinations of as using quenched and 
dynamical Quarkonium simulations. These in- 
clude investigations of the TOg-dependence, deter- 
mination of ttp^""^' from systems and/or split- 
tings with different characteristic scales and in- 
vestigation into the sources of systematic error 
of the ajjg[Mzo] calculation. These tests uphold 

previous determinations of ap and a-jfg and give 
us further confidence in the reliability of those 
results. Our final value for aj^[Mzo] remains 
0.115(2) if the coefficient in the matching to 
MS is set to zero, and 0.117 if this coefficient 
takes its nj = value of 0.96. 
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the Ohio Supercomputer Center and the Atlas 
Centre. We thank the HEMCGC and MILC col- 
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rations and UKQCD, Kogut et al. and Kilcup et 
al. for providing the quenched configurations for 
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